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ABSTRACT 

In  the  development  of  mathematical  software,  often  the  formula  that 
defines  the  mathematical  purpose  of  the  software  is  not  used  directly  in 
the  software.  The  computational  algorithm  used  is  often  mathematically 
equivalent  to  the  defining  formula,  but  bears  little  resemblence  to  it 
for  computational  reasons.  Therefore,  although  the  flow  of  control  of 
the  coded  algorithm  may  be  visible  to  a  maintenance  programmer,  compre¬ 
hending  and  maintaining  the  code  effectively  may  still  be  difficult  if 
only  the  formula  that  defines  the  purpose  of  the  code  is  provided  as 
documentation.  In  this  memorandum,  we  provide  some  examples  of  this 
consequence  of  transforming  the  mathematical  definition  of  a  computation 
into  a  coded  algorithm. 
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INTRODUCTION 

It  is  generally  agreed  that  clearly  written  code  helps  the  pro¬ 
grammer  working  on  it  to  grasp  its  meaning  more  quickly,  so  that  pro¬ 
gram  changes  can  be  applied  with  more  confidence.  Although  structured 
coding  promotes  program  readability  by  exposing  a  program's  flow  of 
control,  so  that  one  can  follow  the  control  flow  in  a  top-down  manner, 
it  does  not  necessarily  follow  that  a  reader  of  structured  code  will 
comprehend  its  algorithms  sufficently  to  effectively  maintain  the  code. 
The  reason  for  this  is  that  in  general  there  are  a  multitude  of  algo¬ 
rithms  that  can  be  chosen  to  perform  any  particular  program  transaction. 
The  algorithm  that  is  chosen  is  determined  by  computational  considera¬ 
tions,  and  often  bears  little  resemblence  to  the  formulas  that  were  used 
to  define  the  transaction  originally.  Therefore,  although  the  coded 
algorithm's  flow  of  control  may  be  highly  visible  to  a  maintenance  pro¬ 
grammer,  comprehending  the  algorithm  may  be  difficult  without  documenta¬ 
tion  that  describes  the  algorithm  and  its  computational  advantages,  as 
well  as  identifying  the  transaction  it  represents. 

In  this  memorandum,  we  give  some  examples  of  mathematically  -  but 
not  computationally  -  equivalent  formulas  that  result  from  transforming 
the  mathematical  definition  of  a  computation  into  a  coded  algorithm, 
showing  that  without  documenting  the  algorithm  the  resulting  struc¬ 
tured  code  may  not  be  comprehensible  enough  to  maintain  effectively. 

BACKGROUND  INFORMATION 

In  the  development  of  mathematical  software,  often  the  formula 
that  defines  the  mathematical  purpose  of  the  software  is  not 
used  directly  in  the  software.  Often  the  computational  algorithm 
is  based  on  some  mathematically  equivalent  formula  that  is  determined 
by  computer  arithmetic,  operating  system  or  hardware  features  that 
impact  computational  accuracy,  execution  efficiency  or  storage  economy. 
For  example,  since  the  laws  of  additive  associativity  and  closure  for 
the  real  number  system  do  not  hold  for  floating-point  computer  number 
systems,  the  mathematically  equivalent  expressions  (x+l)-l  and  x+(l-l) 
will  not  give  the  same  computer  results  for  all  values  of  x. 

Furthermore,  the  same  algorithm  can  be  implemented  in  a  computer 
program  in  different  ways.  For  example,  the  structure  of  the  flow  of 
control  of  a  program  module  depends  on  the  programmer.  In  particular, 
given  two  FORTRAN  implementations  of  an  algorithm,  the  flow  of  control 
of  one  of  them  may  be  easy  to  follow,  while  the  flow  of  control  of  the 
other  may  be  as  difficult  to  follow  as  the  entwining  strands  of  sphaget- 
ti  [1],  On  the  other  hand,  implementation  efforts  to  improve  program  ex- 
cution  efficiency,  which  refine  the  algorithm  further,  may  demote  soft¬ 
ware  clarity. 

Therefore,  three  stages  to  the  development  of  mathematical  soft¬ 
ware  can  be  identified.  Firstly,  a  formula  is  specified  that  defines 
the  purpose  of  the  computation.  Secondly,  a  computational  algorithm  is 
selected  from  among  mathematically  equivalent  forms  of  the  defining 
formula  that  differ  in  their  computational  performance.  This  is  done 
to  produce  an  algorithm  that  satisfies  certain  computational  require¬ 
ments  of  accuracy,  execution  efficiency  and/or  storage  efficiency. 
Thirdly,  the  selected  algorithm  is  coded.  Therefore,  unless  the  map¬ 
ping  of  the  defining  formula  into  the  coded  algorithm  is  documented, 
maintenance  of  the  software  may  still  be  difficult  regardless  of  the 
structuredness  of  the  control  flow  of  the  code. 
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EXAMPLES 

Consider  constructing  software  to  compute  the  magnitude  of  a  com¬ 
plex  number  z  =  x  +  iy  on  a  computer  where  the  largest  and  smallest  pos¬ 
itive  computer  numbers  are  n  and  1/n,  respectively,  with  n  >>  2.  The 
standard  mathematical  definition  for  this  computation  is 

2  2 

abs(z)  =  sqrt(x  +  y  )  (1) 

If  either  the  magnitude  of  x  or  the  magnitude  of  y  is  outside  the  inter¬ 
val  [l/sqrt(n),  sqrt(n)],  use  of  the  standard  definition  as  a  computa¬ 
tional  formula  results  in  computer  overflow  or  underflow.  On  the  other 
hand,  if  the  magnitude  of  x  and  the  magnitude  of  y  lie  in  the  interval 
[1/n,  n/sqrt(2)],  which  contains  the  interval  [l/sqrt(n),  sqrt (n) ], using 
the  mathematically  equivalent  formula 


2 

abs(z)  =  v  sqrt(l  +  (w/v)  ) 


(2) 


where 


v  =  max(abs(x),  abs(y)),  w  =  min(abs(x),  abs(y)) 

does  not  cause  computer  overflow  and  makes  computer  underflow  inconse¬ 
quential.  However,  note  that  in  Figure  1  the  structured  FORTRAN  codes 
based  on  these  mathematically  equivalent  formulas  only  vaguely  resemble 
each  other.  In  fact,  given  that  the  purpose  of  Code  2  is  to  compute  the 
magnitude  of  a  complex  number  z,  which  is  traditionally  defined  by 
Eq.  (1),  it  is  likely  that  without  additional  information  a  maintenance 
programmer  would  have  difficulty  comprehending  the  encoded  algorithm  by 
just  reading  the  FORTRAN  code. 


FIGURE  1 


CODE  1  for  Eq. (1) 


CODE  2  for  Eq. (2) 


ABSZ  =  -1.0 
ROOT  =  SQRT (N/2 . 0 ) 

B  =  ABS ( X ) . LT . ROOT  .AND.  ABS ( Y ) . LT . ROOT 
IF(B)THEN 

ABSZ  =  SQRT (X** 2  +  Y**2) 

END  IF 


ABSZ  =  -1.0 
ROOT  =  N/SQRT (2.0) 

W  =  AMIN1(ABS(X) ,ABS(Y) ) 

V  =  AMAX1(ABS(X) ,ABS(Y) ) 

I F ( V . EQ .0.0) THEN 
ABSZ  =0.0 
ELSE 

I F ( V . LT . ROOT ) 

+  ABSZ  =  V*SQRT ( 1 . + (W/V ) **2 ) 

END  IF 


As  a  second  example,  consider  computing  the  matrix  product  V  of 
the  n-th  order  matrices  A  and  X,  which  is  defined  by 


(V  )  = 
ik 


(T'' A  X 

JTf  «  f 


) 


(3) 
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The  mathematical  definition  of  the  element  in  the  i-th  row  and  k-th 
column  of  V,  given  by  (3),  is  the  inner  product  of  the  i-th 
row  of  matrix  A  with  the  k-th  column  of  matrix  X.  Typically,  one 
computes  all  the  elements  in  the  k-th  column  in  this  way  for  each  of 
the  columns  of  the  product  V.  Expressing  this  in  structured  FORTRAN 
gives 


CODE  3 

DO  30  K  -  1,  N 
DO  25  I  =  1,N 
V(I,K)  =  0.0 
DO  15  J  =  1,  N 

V( I , K )  =  V( I , K )  +  A( I , J)  *  X( J , K ) 
15  CONTINUE 

25  CONTINUE 
30  CONTINUE 


Note  that  for  each  value  of  k  this  algorithm  visits  the  elements 
of  matrix  A  in  row  major  order: 


A  ,  A  , 
11  12 


A  ,  •  •  •  , 
2n 


A  ,  A  , 
nl  n2 


A 

nn 


But,  in  [2],  it  is  shown  that  on  virtual  memory  systems  like  the  VAX, 
addressing  matrix  elements  in  this  order  is  inefficient  when  assigned 
main  memory  is  too  small  to  contain  the  code  and  its  matrices.  In  order 
to  reduce  execution  time,  matrix  algorithms  written  in  FORTRAN  should 
be  based  on  addressing  matrix  elements  in  column  major  order: 


A  , 
11 


A 


f  •  •  •  f 


21 


A  ,  A  , 
nl  12 


A  f  m  m  m  / 
22 


A  ,  .  .  .  , 
n2 


A  ,  •  •  •  , 
2n 


A 

nn 


Now  consider  the  following  mathematically  equivalent  method  of  com¬ 
puting  the  elements  in  the  k-th  column  of  matrix  V,  which  is 
based  on  visiting  the  elements  of  the  A  array  in  column  major  order. 

The  k-th  column  of  V  is  computed  recursively  by  generating  a  finite 
sequence  of  n  +  1  vector  approximations  to  it: 


(0) 

V  =  0, 
k 

where 


(j-1) 

V  +  X  A 

k  jk  j 


( j  =  1,  . . . ,  n)  (4) 


[  A  ] 

C  lj  ] 
[  ] 
[  A  ] 
[  2  j  ] 

[  ] 

A  =  [  .  ] 

j  [  •  ] 

[  .  ] 

[  ] 
C  A  ] 
[  nj  ] 
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(n) 

and  V  is  the  k-th  column  of  V;  hence  at  the  j-th  stage  of  the  re- 
k 

cursion  each  element  of  the  j-th  column  of  matrix  A  is  multiplied  by 
the  j-th  element  of  column  k  of  matrix  X,  so  that  the  algorithm  visits 
the  elements  of  matrix  A  in  column  major  order.  The  corresponding 
structured  FORTRAN  code  for  this  is 

CODE  4 


DO  10  K  =  1,  N 
DO  5  I  =  1,  N 
V(I ,K)  =0.0 
5  CONTINUE 
10  CONTINUE 

DO  30  K  =  1,  N 
DO  25  J  =  1,  N 
DO  15  I  =  1,  N 

V{I,K)  =  V(I,K)  +  A(I,J)  *  X(J,K) 

15  CONTINUE 

25  CONTINUE 
30  CONTINUE 

In  Table  1,  VAX  11/780  batch  execution  times  for  Codes  3  and  4  show 
the  superiority  in  execution  efficiency  of  Code  4  when  matrix 
memory  requirements  exceed  a  users  assigned  main  memory  extent  of  250 
pages.  It  is  assumed  that  V(I,K)  is  accumulated  in  double  precision. 


TABLE  1 


TIMING  FOR  MATRIX  MULTIPLICATION  ON  THE  VAX  11/780 
(assigned  main  memory  extent  =  32,000  words) 


Order 

Memory  Requirements 

Code  3 

Cod 

(n) 

( in  words) 

(  in  cpu  sec. 

)  (  in  cpu 

16 

1,024 

.07 

.07 

32 

4,096 

.61 

.53 

64 

16,384 

5.75 

4.32 

128 

65,536 

86.41 

45.91 

200 

160,000 

910.88 

157.78 

240 

230,400 

5581.4 

274.72 

256 

262,144 

6828.3 

362.44 

512 

1,048,576 

56858.0 

2668.27 

Unfortunately,  knowing  only  the  usual  definition  (3)  of  matrix 
multiplication  and  that  Code  4  computes  the  product  of  two  n-th  order 
matrices,  a  maintenance  programmer  might  replace  Code  4  with 
the  shorter  Code  3  in  order  to  reduce  program  control  complexity. 
However,  doing  this  could  diminish  the  code's  execution  efficiency 
dramatically.  In  other  words,  failure  to  provide  documentation 
that  describes  the  computational  ramifications  of  the  algorithm  used  to 
implement  a  program  module's  function  can  lead  to  counter-productive 
code  modifications  during  software  maintenance. 
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Now  consider  unrolling  the  innermost  DO-loops  of  Code  4  to  reduce 
the  loop  overhead  of  incrementing  the  value  of  I,  testing  the  new  value 
of  I  against  N  and  branching  to  the  beginning  of  the  loop.  Doing  this 
will  improve  the  program's  execution  efficiency  and  refine  algorithm  (4) 
further.  Assuming  that  4  divides  N  exactly,  unrolling  Code  4  to 
a  depth  of  4  gives  the  more  efficient,  but  longer  structured  code 


CODE  5 

DO  10  K  =  1,  N 

DO  5  I  =  1,  N,  4 
V( I , K )  =  0.0 
V( 1+1 ,K)  =  0.0 
V( 1+2 ,K)  =  0.0 
V( 1+3 ,K)  =  0.0 
5  CONTINUE 

10  CONTINUE 

DO  30  K  =  1,  N 
DO  25  J  -  1,  N 

DO  15  I  =  1,  N,  4 

V(I,K)  =  V( I , K )  +  A( I , J )  *  X(J,K) 

V( 1+1 , K )  =  V ( 1+1 ,K)  +  A(I+1,J)  *  X(J,K) 

V ( 1+2 , K)  =  V( 1+2 ,K)  +  A( 1+2 , J)  *  X(J,K) 

V( 1+3 ,K)  =  V( 1+3 ,K)  +  A( 1+3 , J)  *  X(J,K) 

15  CONTINUE 

25  CONTINUE 

30  CONTINUE 


The  corresponding  algorithm  for  computing  the  k-th  column  of  V  is 


(j,i)  ( j -1 , i )  (i) 

V  =  V  +  A  X  (j  -  1, 

k  k  j  jk 

where  for  each  value  of  j 


( j  ,  i) 

V 

k 


[ 

[ 

[ 

[ 

[ 

[ 

[ 

[ 

[ 

[ 


(j) 

V 

4i-3,k 


(j) 

V 

4  i ,  k 


] 

] 

] 

] 

] 

] 

] 

] 

] 

] 


.  .  .  ,  n)  (5) 


[  A  ] 

[  4 i-3 , j  ] 

[  .  ] 
=  [  .  ] 
[  .  ] 
[  A  ] 

[  4  i  ,  j  ] 


(  i  =  l , . . . , n/4 ) 


Note  that  without  an  explanation  of  loop  unrolling,  a  maintenance  pro¬ 
grammer  may  puzzle  over  the  structured  FORTRAN  implementation  of  algo¬ 
rithm  (4)  given  by  Code  5. 


As  a  third  example,  consider  computing  the  eigenvalues  of  a  real 
symmetric  Toeplitz  matrix.  The  elements  of  a  symmetric  Toeplitz 
matrix  T  of  order  n  satisfy  the  relationship 


T  =  T  =  T 

jk  j+l,k+l  kj 


(6) 
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Therefore , 


T  =  T  =  T  (7) 

jk  n-k+l,n-j+l  n-j+l,n-k+l 

by  adding  n+l-j-k  to  both  j  and  k.  Hence  reversing  the  order  of  the 
columns  and  rows  of  a  symmetric  Toeplitz  matrix  yields  the  original 
matrix;  .  in  matrix  notation 

JTJ  =  T  (8) 


where  J  is  obtained  by  reversing  the  columns  of  the  n-th  order  identity 
matrix  I.  A  consequence  of  (8)  is  that  any  real  symmetric  Toeplitz 
matrix  T  of  even  order  2n  can  be  written  in  the  form 


t  t 

where  A  =  A,  B  =  B,  A 
matrix  T  is  similar  to 


[A  BJ  ] 

T  =  [  ] 

[JB  JAJ  ] 

is  Toeplitz  and  B  is  Hankel. 


(9) 


But  note  that 


t  [  A+B  0  ] 

R  TR  =  [  ] 

[  0  A-B  3 


(10) 


where  R  is  the  orthogonal  matrix 

[I  I  3 

R  =  sqrt( .5) [  ]  (11) 

[  J  -J  3 

t 

Therefore,  the  eigenvalues. of  matrices  T  and  R  TR  are  the  same. 
Consequently,  computing  the  eigenvalues  of  an  even  order  real  symmetric 
Toeplitz  matrix  T  is  mathematically  equivalent  to  computing  the  eigen¬ 
values  of  the  smaller  symmetric  matrices  A+B  and  A-B,  where  the  elements 
in  the  k-th  column  of  matrix  A+B  are  given  by 


for  k=l , . . 


(A+B)  =  T  +  T 


jk  j , 2n-k+l 

+  T 

l,k-j+l  l,2n+2-j-k 

if  j 

<  k 

(12) 

+  T 

j-k+1,1  2n+2-j-k, 1 

if  j 

>  k 

Now  consider  constructing  software  that  computes  the  eigenvalues 
of  a  real  symmetric  Toeplitz  matrix  T  of  even  order  n  by  using  a  canned 
subroutine  that  computes  the  eigenvalues  of  any  real  symmetric  matrix, 
where  by  definition  (6) 
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T  =  T  if  j  <  k  (13) 

jk  l,k-j+l 

=  T  if  j  >  k 

If  j-k+1 

Given  the  first  row  TR1  of  matrix  T,  two  possible  structured  FORTRAN 
subroutines  that  compute  matrix  T's  eigenvalues  EIG  using  a  library  sub¬ 
routine  EIGRS  for  the  eigenvalues  of  a  real  symmetric  matrix  are 

CODE  6 

SUBROUTINE  TOPE IG ( TR1 , N , T , EIG ) 

DIMENSION  T(N,N) ,EIG(N) ,TR1(N) 

DO  20  K  =  1 , N 
DO  10  J  =  1 , N 

I F ( J . LE . K ) THEN 
LI  =  K  -  J  +  1 
ELSE 

LI  =  J  -  K  +  1 
END  IF 

T(J,K)  =  TRl(Ll) 

10  CONTINUE 
20  CONTINUE 

CALL  EIGRS  (T,  N,  EIG) 

RETURN 
END 

CODE  7 

SUBROUT I NE  TOPE I G ( TR1 , N , ND I V  2 , APLUSB , E I G ) 

INPUT:  TR1  =  FIRST  ROW  OF  TOEPLITZ  MATRIX 
N  =  ORDER  OF  TOEPLITZ  MATRIX 

NDIV2  =  N/2 ,  THE  ORDER  OF  SCRATCH  MATRIX  APLUSB 
APLUSB=  SCRATCH  MATRIX  OF  ORDER  NDIV2 
OUTPUT:  EIG  =  ARRAY  CONTAINING  THE  N  EIGENVALUES 

DIMENSION  APLUSB ( NDIV2 ,NDIV2 ) ,EIG(N) ,TR1(N) 

SIGN  =  1.0 
DO  40  I  =  1,2 

I POINT  =  (I  -  1 ) *ND I V2  +  1 
DO  20  K  =  1 , NDI V2 
DO  10  J  =  1 , NDI V2 

L2=N+2-K-J 
I F ( J . LE . K ) THEN 
Ll  =  K  -  J  +  1 
ELSE 

LI  =  J  -  K  +  1 
END  IF 

APLUSB (J , K)  =  TRl(Ll)  +  TR1(L2)*SIGN 
10  CONTINUE 

20  CONTINUE 

CALL  EIGRS  (APLUSB,  NDIV2 ,  EIG(IPOINT)) 

SIGN  =  -1.0 
40  CONTINUE 
RETURN 
END 


C 

C 

c 

c 

c 

c 

c 
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Code  6  computes  the  eigenvalues  of  T  directly,  while  Code  7  performs  the 
mathematically  equivalent  computation  of  finding  the  eigenvalues  of  the 
smaller  symmetric  matrices  A+B  and  A-B.  Code  7  is  computationally 
superior  to  Code  6,  since  it  requires  75%  less  storage  for  matrices  and 
executes  75%  faster  when  matrix  T  is  of  sufficiently  large  order  [3,4]. 
However,  knowing  only  definition  (13)  of  a  real  symmetric  Toeplitz 
matrix  T  and  that  Code  7  computes  the  eigenvalues  EIG  of  matrix  T  (if  it 
is  of  even  order  n)  by  using  a  canned  eigenvalue  routine  EIGRS  for  any 
real  symmetric  matrix,  a  maintenance  programmer  may  be  puzzled  by 
Code  7  and  tempted  to  replace  it  with  the  shorter  and  more  comprehen¬ 
sible  (but  less  efficient)  Code  6. 


As  a  final  example,  consider  approximating  the  scaled  Airy  func¬ 
tion  exp(zta)  Ai(z)  for  values  of  z  of  large  magnitude  by  evaluating  a 
partial  sum  of  the  asymptotic  series  [5] 


where 


exp(zta)  Ai (z)  ~ 


-0.5 

( * 


-0.25 
z  /2) 


£ 

Lr  =  o 


k  -k 
(-1)  c  zta 
k 


1.5 

zta  =  z  /1.5,  abs(  arg  z  )  <  n 


-1 

c  =  1,  c  / c  —  k/2  +  5 ( k+1 )  /72 

0  k+1  k 


(14) 


(15) 


Writing 


k  -k 
G  =  (-1)  c  zta 
k  k 


we  see  that  the  terms  of  the  sum  in  (14)  can  be  generated  recursively: 

-1 

G  =  -(c  /c  )  zta  G  (k  =  0,1,2,  ...)  (16) 

k+1  k+1  k  k 

and,  a  fortiori. 


-k-1  k  -k  -k-1 

abs(G  )  =  abs(zta  )  n  (c  /c  )  =  2  k!  abs(zta  )5/72  (17) 

k+1  j=0  j+1  j 


For  any  specific  value  of  z,  the  approximation  to  the  scaled  Airy  func¬ 
tion  having  maximum  accuracy  is  obtained  by  evaluating  the  partial  sum 


(18) 


where  G 

n+1 


is  the  first  term  encountered  in  the  series  for  which 

abs(G  )  >  abs(G  ) 
n+1  n 
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By  (16),  condition  (19)  is  equivalent  to 

c  /c  >  m  =  abs(zta)  (20) 

n+1  n 

Therefore,  the  first  value  of  n  for  which  (19)  holds  is 

2  0.5  -1 

n  =  [(9m  +  9m  +  1)  /3  +  m  -  0.5]  ~  [2m  -  (14.4m)  ]  (21) 

where  [x]  is  the  smallest  integer  greater  than  or  equal  to  x.  How¬ 
ever,  if  m  is  sufficiently  large  computation  of  G  will  cause  com¬ 
ic 

puter  underflow  well  before  k  =  n,  and  the  evaluation  of  (18)  out  to 

G  then  becomes  impractical.  Consequently,  in  a  practical  computation, 
n 

(18)  should  be  evaluated  out  to  G  or  the  G  with  the  smallest  magni- 

n  k 

tude  that  does  not  cause  computer  underflow,  whichever  comes  first. 


The  largest  value  of  m  for  which  (18)  can  be  summed  out  to  G  with- 

n 

out  any  of  the  G  underflowing  depends  on  the  smallest  positive  computer 

k 

number  s.  A  good  approximation  to  this  value  of  m  is  obtained  by  set¬ 
ting  the  approximation  (17)  for  abs(G  )  equal  to  s  with  k  =  2m,  ap- 

k+1 

plying  Stirling's  asymptotic  formula  for  factorials  and  then  solving  the 
resulting  equation  for  m.  This  gives  the  equation 

-0.5 

m  =  -ln(m)/4  -  ln(7.2sTr  )/2  (22) 

whose  solution  is  m* ,  the  limit  of  the  rapidly  convergent  sequence 

-0.5 

m  =  -ln(m  )/4  +  m  ,  m  -ln(7.2s7r  )/2  (23) 

j  j-1  0  0 

In  other  words,  for  any  given  value  of  z  for  which  m  does  not  exceed 

m*,  (18)  can  be  summed  out  to  G  without  any  of  the  G  underf lowing ; 

n  k 

if  m  exceeds  m* ,  (18)  can  be  summed  out  to  the  term  whose  index  k  is 
-1 

[2m*f  -  (14.4m)  ]  just  short  of  G  underflowing  for  a  positive  value 

k 


of  f  <  1. 
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To  find  f  we  proceed  in  the  same  way  that  we  did  to  find  m*,  but 

we  set  the  approximation  (17)  equal  to  s  with  k  =  2m*f.  This  gives 

f  =  (• (f  +  l/4m*)ln(f)  -  c  )  /  ln(em/m*)  (24) 

c  =  (l/2m*)ln(  7.2sm  /  sqrt(7rm*).  ) 

whose  solution  is  the  limit  of  the  sequence 

f  =  (  (f  +  l/4m*)ln(f  )  -  c  )  /  ln(em/m*)  (25) 

j  j-1  j-1 

for  an  initial  value  f  <  1.  If  f  =  l/log(m),  then  f  gives  an  adequate 

0  0  2 

approximation  to  f. 


To  summarize,  evaluating  the  asymptotic  series  until  (19)  is  sat¬ 
isfied  or  G  underflows  is  similar  to  evaluating  the  series  out  to  the 
k  -1 

term  whose  index  k  has  the  value  [ 2m*f- ( 14 . 4m)  ],  where  f  =  m/m*  if  m  <  m* 

For  any  given  z,  the  latter  method  provides  an  apriori  estimate  of  the 
optimal  number  of  terms  of  the  series  to  sum,  eliminating  the  need  for 
the  comparison  test  (19)  during  the  summation  of  the  series. 

Therefore,  two  possible  structured  FORTRAN  subprograms  for  evalu¬ 
ating  the  asymptotic  series  are 

CODE  8 

FUNCTION  SUM  ( ZTA ) 

COMPLEX* 16  SUM,  ZTA,  GK,  RZTA 

DOUBLE  PRECISION  ABSGK ,  ABSGK1 ,  S,  FACTOR 

DATA  S/2.94D-39/ 

SUM  =  DCMPLX ( 0 . 0D0 ,  0 . 0D0 ) 

K  =  0 

GK  =  DCMPLX( 1 . 0D0 , 0 . 0D0 ) 

RZTA  =  GK  /  ZTA 
ABSGKl  =  1.0D0 

C  DO  UNTIL  (  ABSGKl. GE. ABSGK  .OR.  ABSGKl. LE.S) 

5  SUM  =  SUM  +  GK 

FACTOR  =  0.5D0  *  K  +  5.0D0  /  (K  +  1)  /  72.0D0 

K  =  K  +  1 

GK  =  -  FACTOR  *  RZTA  *  GK 
ABSGK  =  ABSGKl 
ABSGKl  =  CDABS(GK) 

IF  (  ABSGK1.LT. ABSGK  .AND.  ABSGKl . GT . S ) GO  TO  5 
CONTINUE 
RETURN 
END 
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CODE  9 

FUNCTION  SUM  (ZTA) 

COMPLEX* 16  SUM,  ZTA,  GK ,  RZTA 
DOUBLE  PRECISION  FACTOR 

REAL  MSTAR,  EOMSTR ,  TMSTAR ,  R4MSTR,  CLN,  M,  F,-C,  D 
DATA  MSTAR/ 42.72/,  TMSTAR/8 5 . 44/ ,  EOMSTR/ . 063630 19/ 

DATA  CLN/-89. 198027/,  R4MSTR/5 . 8520599 E- 03/ 

M  =  CDABS ( ZTA ) 

I F ( M  .GT.  MSTAR) THEN 
F  =  1 . 0/ALOG10 (M) 

C  =  (CLN  +  ALOG(M) ) /TMSTAR 
D  =  ALOG ( EOMSTR*M ) 

F  =  ( ( F  +  R4MSTR ) *ALOG ( F )  -  C)/D 

F  =  ( ( F  +  R4MSTR ) * ALOG ( F )  -  C)/D 

N  =  2.0  *  (MSTAR  *  F  +  0.5)  -  .069444444  /  M 

ELSE 

N  =  2.0  *  (M  +  0.5)  -  .069444444  /  M 
END  IF 

SUM  =  DCMPLX ( 1 . 0D0  ,  0.0D0) 

GK  =  SUM 

RZTA  =  GK  /  ZTA 

DO  100  K  =  0,  N-l 

FACTOR  =  0.5D0  *  K  +  5.0D0  /  (K  +  1)  /  72.0D0 
GK  =  -  FACTOR  *  RZTA  *  GK 
SUM  =  SUM  +  GK 
100  CONTINUE 
RETURN 
END 

Code  8  terminates  summation  of  the  series  when  either  (19)  is  satis¬ 
fied  or  a  term  of  the  series  underflows  s,  while  Code  9  terminates 
when  an  apriori  estimate  of  the  summation  index  is  reached  for  which 
either  one  of  these  conditions  is  satisfied.  Although  Code  9  is  longer 
than  Code  8,  the  VAX  11/780  executes  Code  9  30%  faster  than  Code  8; 
in  fact,  if  an  initial  value  for  the  sequence  (25)  were  found  so  that 
f  would  give  an  adequate  approximation  to  f,  then  Code  9  would  execute 
1 

approximately  40%  faster  than  Code  8.  In  any  event,  knowing  only  that 
the  purpose  of  these  codes  is  to  evaluate  the  sum  (18)  for  a  given  value 
of  zta  until  either  (19)  is  satisfied  or  G  underflows  s,  without  ad- 

k 

ditional  information  a  maintenance  programmer  would  have  difficulty 
comprehending  how  this  is  accomplished  by  the  more  efficient  Code  9. 


CONCLUSION 

Although  structured  coding  promotes  software  clarity  by  highlight¬ 
ing  the  flow  of  control  of  the  code,  the  result  of  searching  for  a  com¬ 
putational  method  that  hopefully  performs  the  same  function  optimally 
may  demote  software  comprehensibility.  Since  the  algorithm  that  is 
finally  implemented  can  depart  significantly  from  the  formula  that  was 
used  originally  to  define  the  purpose  of  the  code,  the  algorithm  must 
be  documented  sufficiently  or  else  the  code  may  be  difficult  to  main¬ 
tain  effectively. 
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